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CONSISTENT RECALIBRATION OF YIELD CURVE MODELS 


PHILIPP HARMS, DAVID STEFANOVITS, JOSEF TEICHMANN, AND MARIO V. WUTHRICH 


Abstract. The analytical tractability of affine (short rate) models, such as the Vasicek 
and the Cox-Ingersoll-Ross models, has made them a popular choice for modelling the 
dynamics of interest rates. However, in order to account properly for the dynamics of 
real data, these models need to exhibit time-dependent or even stochastic parameters. 
This in turn breaks their tractability, and modelling and simulating becomes an arduous 
task. We introduce a new class of Heath-Jarrow-Morton (HJM) models that both fit the 
dynamics of real market data and remain tractable. We call these models consistent re¬ 
calibration (CRC) models. These CRC models appear as limits of concatenations of for¬ 
ward rate increments, each belonging to a Hull-White extended affine factor model with 
possibly different parameters. That is, we construct HJM models from “tangent” affine 
models. We develop a theory for a continuous path version of such models and discuss 
their numerical implementations within the Vasicek and Cox-Ingersoll-Ross frameworks. 


1. INTRODUCTION 


1.1. Principles of yield curve modelling. Modelling the stochastic evolution of yield 
curves is an important task in risk management, forecasting, decision making, pricing and 
hedging. We emphasise here three principles of yield curve modelling (or any other traded 
instrument in finance): we certainly require that all models for traded assets’ prices are 
free of arbitrage; therefore we do not state this as a principal requirement. 


Robust calibration: the model is selected simultaneously from time series and 
prevailing market prices, as explained in 12 . Model parameters which are invariant 
under equivalent measure changes should be estimated by a statistical procedure 
from time series data. The remaining parameters are calibrated by solving an 
inverse problem with respect to the prevailing market prices. All model parameters 
should be constant during the life time of the model; only state variables may 
change. 

Consistency: an interest rate model is called consistent if the stochastic process 
of yield curves does not leave a pre-specified set I of possible market observables 
(in 16 the set X is assumed to be a finite dimensional sub-manifold of curves 


corresponding to a curve fitting method). Here, we add the following requirement: 
the yield curve process should be able to reach any neighbourhood of any yield curve 
in X with positive probability because any newly arriving market configuration 
is a possible model state. Consequently, the model can be recalibrated to a new 
market configuration without losing consistency to the model choice with old 
parameters; we say that the model satisfies the consistent recalibration property. 
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• Analytic tractability: relevant quantities of a model can be calculated quickly 
and accurately. In particular, one should be able to simulate state variable 
increments efficiently. This can be a delicate problem in the presence of boundary 
conditions. 

We briefly comment on some of these principles. 

• By a model for the term structure of interest rates, we understand a fully specified 
stochastic process taking values in the pre-specified set of yield curves I. We shall 
always consider a parametrised class of models consisting of one fully specified 
model for each initial state in I and each parameter value. 

• In practice, interest rate models are recalibrated on a regular basis (e.g. daily) to 
market data. Suppose that the consistent recalibration property does not hold 
for today’s model. Then tomorrow’s market yield curve might lie outside of the 
set of possible realisations of today’s model. If this happens, then tomorrow’s 
recalibration necessarily implicates a rejection of today’s model. On the other 
hand, no inconsistencies occur if recalibration is an update of the state variables 
of the model and does not involve a change of model parameters. 

• The robust calibration principle separates the easier task of estimating volatilities 
from the more difficult task of estimating drifts. Moreover, it tells us exactly which 
parameters may be estimated from time series (namely those which are invariant 
under equivalent changes of measure). In [22| we use our results to model and 
filter the market price of risk. 


1.2. Consistent recalibration models. The goal of this work is to present a new model 
class for yield curve evolutions satisfying all three principal requirements with respect to 
sets I which are sufficiently large to be of practical use (think of open subsets of a Hilbert 
space of curves). Mathematically speaking, we look for yield curve models with full or 
large support, which are in addition analytically tractable. 

Often the full support property does not accord with analytic tractability beyond 
elliptic models, which are too restrictive in infinite dimension. We illustrate this with 
an example. Consider the Hull-White extended Cox-Ingersoll Ross (CIR) model. In this 
model the short rate is given by the SDE 

dr{t) = {0{t) + fdr{t))dt -\- a/ ar(t)dW(t), 


where 9{t) >0 determines the time-dependent level of mean reversion, /3 < 0 the speed 
of mean reversion, and a > 0 the level of volatility. 

The model can be calibrated to any initial yield curve from a large subset I of curves 
by choosing an appropriate Hull-White extension 9. For any fixed initial yield curve, the 
distribution of yield curves at some future t > 0 is concentrated on a one-dimensional 
affine subspace of curves. Therefore, market observations are generally not in the support 
of the model, and the consistent recalibration property does not hold with respect to 
I. Moreover, the low dimensionality of the model is apparent at the level of realised 
covariations of yields for different times to maturity: the matrix of covariations has rank 
one, which is in stark contrast to observations from the market (see [Figure 7.13 ). Finally, 
calibrated model parameters vary significantly over time as shown in Figures |7.4| and |7.6[ 
which contradicts the requirement of robust calibration. 

As a remedy, one could make some model parameters stochastic and include them as 
state variables. For example, one could make a = ay and (3 = fdy depend on a parameter 
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y and write dynamics of the form 

dr{t) = {0{t) + l3Y{t)r(t))dt + sj aY{t)r{t)dW{t), 
dY{t) = fj.(Y{t)^dt + a(Y{t)^dW{t). 

Unfortunately, this usually breaks the analytic tractability of the model in the sense that 
zero-coupon bond prices cannot be calculated anymore analytically. 

The key idea of consistent recalibration (CRC) models is to lift the short rate model 
to a HJM model and to introduce stochastic parameters on that level. Let h{t) denote 
the forward rate curve at time t in Musiela parametrisation (i.e., as a function of time 
to maturity). Then CRC models are defined by the joint dynamics 

dh{t) = (^Ah{t) + + <7^l^{r{t))dW{t), 

dY(t) = y,{Y(t))dt + a {Y (t)) dW (t), 


where r(t) = h(t)(0), A is the generator of the shift semigroup, and are the 

HJM drift and volatility of the short rate model with constant parameter y. This model, 
and more generally the class of CRC models, has the following properties: 


It is a full-fledged HJM model providing the benefits of robust calibration. In¬ 
deed, all model parameters can be estimated from realised covariations of yields 
rather than calibrated by solving high-dimensional inverse problems. Thus, the 
parameters are estimated from yield curve dynamics instead of calibrated to static 
yield curves, while an exact match to the current yield curve is guaranteed by the 
Hull-White extension 6. 

In the language of term structure equations, see [^[^[^, CRC models are tangent 
to affine factor models. This means that they are (limits of) concatenations of affine 
factor models. The construction is illustrated in [Figure The concatenated 
models are allowed to have distinct static parameters. Making the parameters 
stochastic (in an independent or dependent way) will lead generically towards con¬ 
sistent recalibration, since the conditions for finite-dimensional realisations are not 
fulfilled anymore. The consistent recalibration property is also reflected in the much 
higher ranks of the covariation matrices of yields with different times to maturity. 
Indeed, our empirical analysis shows that they are closer to those observed in the 
market than in the corresponding models without CRC extension (see Figure 7.131. 
Despite all this flexibility, the model remains analytically tractable', zero-coupon 
bond prices are state variables, and state variable increments can be simulated effi¬ 
ciently because they look infinitesimally like Hull-White extended affine processes. 
This means that the Fourier transform of the infinitesimal increments is known 
and that all sampling techniques of affine processes apply. In particular, efficient 
high-order positivity-preserving simulation schemes for the CIR process such as 
can be used. No similar schemes are available for general HJM equations 
with non-Lipschitz vector fields. In our numerical implementation, we achieve 
first-order convergence of the splitting scheme. We also give a theoretical proof 
of first-order convergence in the Vasicek case. 


These properties are important in risk management and in the current regulatory 
framework , where one needs tractable and realistic (non-Gaussian) models of long-term 
returns on bond portfolios. Moreover, the CRC approach allows one to easily implement 






4 PHILIPP HARMS, DAVID STEFANOVITS, JOSEF TEICHMANN, AND MARIO V. WUTHRICH 


stress tests for risk management purposes by selecting a suitable model for the parameter 
process. First evidence of improved fits is provided in 


22 


The principle behind CRC models applies also to the modelling of more general term 
structure dynamics. For example, it could be applied to multi-curve interest rate models 
and to models of the term structure of option prices 


29 


1.3. Organisation of the paper. In Section 2 we discuss relations to other yield curve 
modelling approaches. In [Section 3i we introduce Hull-White extended affine short rate 
models, which are the building blocks of CRC models introduced in lSectioiTd In Sectionsj^ 
and 1^ the one-factor Vasicek and CIR case is developed in full detail. In Section 7] our 
numerical implementation and some empirical results are presented. 

2. Relations to other yield curve models 

Several in part overlapping approaches to yield curve modelling have been developed. 
The models can be roughly categorised as factor, HJM, principal component analysis 
(PCA), and filtered historical simulation models. We briefly analyse these models with 
respect to our requirements and compare them to the new class of CRC models. 

2.1. Factor models. Factor models are based on a factor process, which usually describes 
certain market factors, from which - by means of basic principles - the entire yield curve 
can be derived (see 17 for an overview). Let X = (X(t))t>o be a factor Markov process 


acting on a finite dimensional state space and depending on a parameter vector y, and let 
B := g(X) be the bank account process, for some positive functional g. Then one obtains 
- with respect to the pricing measure P - the relation 

B(t) 


P{t,T) =E 


B{T) 


Ht) 


= G(t,r,A(t)), 


for some function G also depending on the parameter vector y. 

Market data arrive in the form of daily yield curves. By means of calibration the 
initial state xq and a parameter vector j/o are chosen to explain today’s market data. 
By choosing the parameter vector rich enough, one receives good fits to today’s market 
data. Apparently a recalibration at time t = \ can (and will) lead to another state xi and 
another parameter vector yi. As states may vary stochastically, the change of xq to Xi 
is in principle not a problem, but the change of parameter vector is. This means that one 
has to decide at time t = 1 whether to continue with the model specified by j/o whether 
to switch to the model specified by yi. This problem can be alleviated to some extent 
by using a combination of filtering and calibration techniques to stabilise the choice of y, 
as described in, e.g., 19 Nevertheless, robust calibration remains an unresolved issue. 


The consistent recalibration property does not hold unless the set I is very small. (If I 
is a sub-manifold, its dimension cannot be larger than that of the factor process.) However, 
on the positive side, factor models are often analytically tractable, for instance, within 


the affine class (see e.g. 17 15]). 


2.2. HJM-models. Markovian HJM models are an extreme version of factor models: the 
yield curve itself is taken as state variable (possibly together with some hidden state vari¬ 
ables). Calibration to daily arriving yield curves is now a matter of statistical estimation 
from the time series of market data. An appropriate parametrisation of instantaneous 
co-variance, jump structure, and drifts will lead to a statistical inference problem, an 
infinite dimensional one though. Hence, the paradigm of robust calibration, including 
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the requirement of calibration through estimation, is fulfilled in the optimal sense. If the 
process acting on yield curves is “irreducible”, i.e., every neighbourhood of a state can be 
reached with positive probability, then even consistent recalibration is possible (see [^). 
However, one usually encounters a severe lack of analytic tractability within this model 
class. Euler and higher order schemes (often) require strong assumptions on the vector 
fields (c.f. [^). Usually no exact or high-order simulations of infinitesimal increments 
are at hand in contrast to CRC models, where this is often the case. 


2.3. PCA- or local PCA models. Principal component analysis (PCA) or local PCA 
considers yield curves as outcomes of a statistical model, which is estimated by standard 
PCA techniques (see e.g. [^[^[^). When the statistical model is too simplistic, often 
arbitrage enters the field, which is an undesirable feature. A more refined version is actually 
equivalent to a HJM model with constant vector fields (as e.g. in |^). Here preserving 
floored interest rates, which is desired in some situations, is not possible. PCA inspired 
models, correctly implemented, allow for robust calibration and consistent recalibration, 
but are usually not very tractable from an analytic point of view. 


2.4. Filtered (historical) simulation. Historical simulation is a standard industry 
technique to simulate distributions of yield curves by considering the relative returns as 
independent samples of an unknown distribution, see 19, ^ . Certainly this assumption 


can lead to difficulties with the absence of arbitrage, but this can be solved as in 28 32 


The most important problem is the state-independence of the distribution. Again also 
(filtered) historical simulation can be embedded into the realm of HJM models. These 
models then allow for robust calibration and consistent recalibration, but are usually not 
very tractable from an analytic point of view. 


3. Hull-White extended affine short rate models 


3.1. Overview. We set the stage for CRC models by describing Hull-White extended affine 
short rate models, focusing first on the correspondence between Hull-White extensions and 
initial forward rate curves. The one-dimensional short rate model of the introduction is 
replaced by more general multi-dimensional factor models for the short rate. The parameter 
y, which becomes stochastic in the CRC setting, is kept constant and fixed for the moment. 


3.2. Setup and notation. (H, (J^(t))t>o, P) is a filtered probability space. The filtra¬ 

tion satisfies the usual conditions. The measure P plays the role of a risk-neutral measure. 
All processes are defined on U, adapted to (J^(t))t>o, and cadlag. W = (lU(t))t>o is a 
d-dimensional (J^(t))t>o-Brownian motion. 

The short rate process r = {r{t))t>o is determined by a factor process X = {X{t))t>o 
with values in a state space X. The evolution of the factor process depends on a parameter 
process Y = {Y{t))t>o with values in a space Y. In all of Section 3 the parameter process 


Yit) = y is assumed to be constant and fixed, whereas it is allowed to vary in Section 4 
below. 

The spaces X and Y are both subsets of some finite dimensional vector space. X is, up 
to permutation of coordinates, of the canonical form x with di -I- d 2 = d> 1. The 
canonical basis vectors in are denoted by Ci,..., Cd, and (•, •) denotes the Euclidean 
scalar product. Of course we could consider more general affine processes here, which 
take values, e.g., in products of cones of positive-semidefinite matrices and real lines like 
Wishart-Heston models. 
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For each (x, y) S X x Y, there is a symmetric positive semidehnite matrix Ay{x) G 
and a vector By{x) G determining the volatility and the drift of X. The expressions 
Ay{x) and By(x) are affine in x, i.e., 




By(x) = by+j2^y 


for all (x, y) € X X Y, 


i=l 


i=l 


for symmetric positive semidehnite matrices Uy , a^, ... ,ay G and by, Py, ..., jSy G 
We denote by ^yAy(x) the symmetric positive semidehnite square root of Ay(x). Note: 
other choices of square roots are possible, and [Assumption 3.1 below does not depend on 
the choice of square root by 30 Chapter V. 19-20]. Moreover, a function 0 G C(R+) is 


given, which is used to make the drift of X time-inhomogeneous. 


3.3. Factor process and short rate. The factor process A is a continuous, X-valued 
solution of the SDE 

(3.1) dX(t) = Ay{X{t))dW(t) + {0{t)ei + By{X{t))^dt 
with initial condition A(0) = x S X. The short rate is given by 

(3.2) r{t) = {\,X{t)), for alH > 0, 

for some hxed £ gR and X gR‘^ satisfying (A, ei) 7 ^ 0. 


Assumption 3.1. It is assumed that SDE (3.1) has a unique continuous 
X, for each initial condition A(s) = x, where (s, x) 


-valued solution 
+ X X. In this case, the parameters 


(y, 9) are called admissible. Moreover, it is assumed that X satishes the moment condition 


(3.3) 


E 


,-/o(^+(AAC)»ds 


< 00, 


for all t > 0 . 


3.4. Exponential moments and Riccati equations. The process X, or rather the 
family of processes obtained by varying the initial conditions in SDE (3.1), is time- 
inhomogeneous affine. All coefficients in SDE (3.1) are independent of time, except for 
the drift 0; we call X Hull-White extended affine and 9 its Hull-White extension. This we 
are going to highlight in detail below. Our main reference for time-inhomogeneous affine 
processes is [T^ . 

Functions ($y,'I'y) G C°°(M+)xC°°(IR+;K'^) are called solutions of the Riccati equations 
(with parameter y gY) if 


(3.4a) % = Fy° ‘J’y(O) = 0, 

(3.4b) % = RyO^>y- A, ^/.(O) = 0 

holds, where {Fy,Ry) G C{R‘^) x K"^) are given by 

Pyiu) = ^{u,ayu) -G {u,by), Rfyffi) = ]^{u,alu) + {u,Pl), 

for all u gR'^ and i G {1,..., d}. 


Lemma 3.2. Let X be a solution of SDE (3.1) for some admissible parameters {y,9). 
Then X satisfies moment condition (3.3) if and only if there exists a solution ($^,4'^) of 
the Riccati equations (3.4) for all times t > 0. Moreover, if there exists a solution, even 
a local one, of the Riccati equations, it is unique. 
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Note that the Riccati equations (3.4) only depend on y, and not on the choice of the 
Hull-White extension 6. 

Proof. Let Z be the X^-valued process Z{t) = {Zi{t), Z 2 {t)) = (X(t), f* X(s)ds). Then 
Z is an Ito diffusion whose drift and volatility at time t > 0 are given by 


(0(t)ei + By(Zi(t)), Zi(t)j € 


r)2d 


Ay{Z,it)) 0 


0 


0 


T>2dx2d 


respectively. Clearly, these expressions are affine in Z(t). Moreover, the time-inhomogeneity 
0{t)ei is (by definition) continuous in t. Therefore, the process Z, or rather the family of pro¬ 
cesses obtained by varying the initial condition of Z, is strongly regular affine, see 18 Theo¬ 
rem 2.14]. For each (t, ui, U 2 ) € K+ the functional characteristics of Z are given by 

F{t,ui,U 2 ) = 9{t){ui,ei) P Fy{ui) G K, 

R{t,ui,U2) = {Ry{ui) + U2,0) G X 


Moment condition (|3.3|), expressed in terms of Z, reads as follows: 

„-(A.Z2(T))' 


E 


< 00, 


for all T > 0. 


By 21 , the moment condition is equivalent to the existence of a solution 1 , 1 ^ 2 ) of 
the following Riccati system associated to Z: 


T) = e{t) (V^i if, T), ei) + Fyd;, (t, T)), 

= Ry('ifi{t,T)) + ijj2{t,T), 

-dtf^2{t,T)=0, 

Equivalently, the relations 'ip 2 {t,T) = —X and 

T) = 0{s){'ify{T - s), ei)ds + %{T - t), 


(/)(T,r) = 0, 
i/-i(T,r) = o, 
MT,T) = -X. 


Mt,T) = ^y{T-t) 


hold identically, where (<l>j^,ff'y) is a solution of the Riccati equations (3.4). Uniqueness 
holds for these equations because the vector fields are locally Lipschitz. □ 

3.5. Bond prices and forward rates. By no-arbitrage arguments zero-coupon bond 
prices in the short rate model of [Section 3.3|are given by 


P(t,T) = E 


r{s)ds 


m 


= E 


Sf (t+XX{s)))ds 


Fit) 


T>t>0. 


For essentials on short rate models we refer to 17 Chapter 5]. We define the (instantaneous) 
forward rates by 

h{t,T) = h(t)(T) = -dr log (p(t,t -I- r)), t, r > 0. 

The parametrisation of the forward rate as a function of t and r is called Musiela 
parametrisation. It is particularly useful in this paper since (h(t))(>Q will be interpreted 
as a stochastic process taking values in a suitable space of functions on K+. 

By the affine nature of the factor process X, zero-coupon bond prices and forward rates 
can be obtained by solving the Riccati system of ODE’s (3.4). 



















8 PHILIPP HARMS, DAVID STEFANOVITS, JOSEF TEICHMANN, AND MARIO V. WUTHRICH 


Theorem 3.3 (Zero-coupon bond price and forward rate). Let X satisfy moment condi¬ 
tion (3.3) and let {^y, ipy) be the unique solution of the Riccati equations with parameter 
y, given by Lemma 3.^ Then the bond prices in the short rate model (3.11-(3.2) satisfy 

(3.5) 

log(P(t, T)) = -£{T -t) + &is){^y{T - s), ei)ds + $,,(T - t) + {^y^T - t), X(t)), 

for all T > t > 0, and the forward rates are given by 

(3.6) h{t,T)=i- [ 0{t + s){%{T-s),ei)ds-<£>'yiT)-{^'y{T),X(t)), 

Jo 

for all t,T > 0. 


Proof. We borrow from the proof of Lemma 3.2 where moment condition (3.3) was shown 
to be equivalent to the existence of solutions of the Riccati system associated 

to the process Z = (X,J X). Moreover, {cf, if 1 , 1 ^ 2 ) are closely related to the solutions 
($y, tby) of Riccati system ( |3.4[ ). By the main theorem in and its conditional version, 
the affine transform formula 


E 


,-(A.Z2(T)> 


Pit) 




for all r > t > 0, 


holds. A direct calculation shows this formula to be equivalent to formula (3.5) for 
bond prices. Formula (3.6) for forward rates is obtained by taking the logarithm and 
differentiating with respect to r. □ 


3.6. Heath-Jarrow-Morton equation. The evolution of forward rate curves is de¬ 
scribed by the HJM equation. For each (a;, y) S Xx Y, let and be given by 

/i™(x) = {%,Ayix)%) e C°°(K+), a™(a;) = -^Ay{P)% G C^{R+-X)- 

Note that the familiar HJM drift condition holds: 

(3.7) Tf^{x){T) = ^a™(x)(r), £ aH™(x)(s)ds^ , for all r > 0. 

Let El be a Hilbert space, destined to contain the forward rate curves of the model. By 
abuse of notation the symbol h is used interchangeably to denote an element of El and 
the forward rate process. 


Assumption 3.4. El is a Hilbert space with the following properties: 

(i) El C C(K+) and the evaluation map eval^: h 1 —)■ /i(r) is continuous on El, for each 
T G M+; 

(ii) for each (x, y, z) G X x Y x fj,^^^ix) and {a^'^^{x), z) are elements of El; 

(in) the right shifts (5(t))(>o mapping h to hit+-) define a strongly continuous semigroup 
on El with infinitesimal generator A. 


Hilbert spaces of forward rate curves which comply with the requirements of |Assump-| 

Sections 5, 7.4.1 and 7.4.2] for the Vasicek and CIR models. 
4) (c.f. “ 


tion 3.4 are constructed in 16 


In the domain I? (A) C El H 


16 


Lemma 4.2.2]) of the infinitesimal generator 


A we can characterise the process (h. A) as follows. 
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Theorem 3.5 (HJM equation). Let {h,X) be given by Theorem 3.3 and assume that 
h(t) S ^{A), for each t>0. Then the process {h,X) is a strong solution of the following 
SPDE on H X X; 


dh{t) = +/i™(X(t)))dt+ cr™(X(t))rftT(t), 

dX{t) = ^ Ay{X{t))dW{t) + (6»(t)ei + 


(3.8) 


Outside of the domain of A the forward rate process can be characterised as a mild 
solution of Equation (3.8) For the concepts of mild, weak, and strong solutions of SPDEs 
we refer to w Section 6.1]. 

In the one-factor case, the factor process {X{t))t>o is a simple functional of the forward 
rate. Then Equation (3.8) can be rewritten as an evolution equation for the forward rate 
process alone (c.f. Equation (5.4)). This is also possible in the multi-factor case, but the 
corresponding functional is more complicated, which is why we choose to present the HJM 
equation in the form (3.8). 


Proof. Differentiating formula (3.6) for forward rates with respect to t and r and using 
'i's/(0) = one obtains for each r > 0 

dh(t,T) = ^ ^ 0(s)(4'"(t -I- r - s), ei)ds + 9(t + T)(A,ei) 

+ 0(f)(<(T), ei) Vt - (<(r), dX(t)), 


/ t + T 

9{s){^y{t + T - s), ei)ds -I- 9{t + r)(A,ei) 
-$"(r)-(iIi"(r),X(t)). 

Subtracting the equations and cancelling out the integral as well as the term next to it yields 
dh{t,T) = {Ah{t,T) + $"(t) -H (4'"(r),X(t)) -H6»(t)(4']^(T),ei)) dt - (^'^(r), dX(t)). 


When dX{t) is replaced by the right-hand side of SDE (3.1), the 0(t)-term cancels out 
and one obtains for each r > 0 

dh{t, r) = {Ah{t, t) + $"(r) + {%{T),Xit)) - {%{T),By{Xm) dt 
- l^^'y{r),^Ay{X{t))dW{t)y 

The symmetric matrix ^ Ay{X[t)) can be moved to the other side of the scalar product, 
and one immediately recognises the volatility a^^^{X{t)). A direct calculation shows 
that the drift is equal to /i^'^'^(A(t)). Indeed, 


+ (4-", x) - (4-;, By{x)) = F;^o^SJy 4-; + { R'y o ^iJy ■ 4-;, x) - (4-;, By{x)) 

= = Aty™(a:), 

which follows from the relations 

Fyiu) ■v = {u, ayv) + {v, by), (A(,)'(u) ■ v = {u, alyv) + {v, P^). 


□ 
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3.7. Forward rates and Hull-White extensions. Relation (3.6) between the forward 
rate h{t) and the Hull-White extension 6 plays a key role in calibration (and recalibration) 
of the model. It can be expressed concisely as 

h{t) ='Hy{S{t)9,X{t)), S{t)9 = Cy(h{t),X{t)), for alH > 0, 

where S{t) is the right shift operator, see Assumption 3.4[pli] ), Hy calculates the initial 
forward rate curve from the Hull-White extension given parameter y, and Cy performs the 
inverse operation of calibrating a Hull-White extension to an initial forward rate curve. 
Formally, for each (t,x,9) S K+ x X x (7(11+), these operators are given by 


S(t)9 = 9(t + -) e C'(K+), 

ny(9, x)=l- Iy{9) - - {%, x) e C'i(M+), 

Iy{9)= [ 9{s){'S'y{--s),e,)dseC\R+). 

Jo 

Note that Hy involves the Volterra integral operator Xy. The operator Cy (the letter C stand¬ 
ing for calibration) is defined as the partial inverse of Uy given by the following theorem. 

Theorem 3.6 (Calibration to initial forward rate curves). Let {h, x) G (7^(11+) x X satisfy 
h{0) = £ + {X,x). Then the Volterra integral equation h = 'Hy{9,x) has a unique solution 
9 G (7(1^+), which we denote by Cy{h,x). 

The theorem is a direct consequence of the following lemma. 

Lemma 3.7. For each y G Y, the Volterra integral operator 

Xy-. (7(K+) ^ {hGC\R+): h{0)=0} 

is bijective. 

Proof. This follows from Theorem 2.1.8], noting that the integral kernel 

(3.9) Ky{s,t) = {'^'yft — s),ei) , for alH > s > 0, 

satisfies \Ky{t,t)\ = |(A, ei)| > 0 and both Ky and dtKy are continuous. □ 

Note that calibration of a Hull-White extension 9 requires the inversion of the Volterra 
integral operator Xy. Here the assumption (A,ei) 7 ^ 0 is needed. 

3.8. Numerical solution of the Volterra equation. In the absence of analytical for¬ 
mulas, the Volterra equation has to be solved numerically. We are aiming at a second 
order approximation to keep the global error of the simulation scheme of order one. Thus, 
we approximate the Volterra integral operator Xy by the trapezoid rule, which yields an 
operator Xy given by 

%{9){t.^) = S l^^('^y(T„),ei}9(0)+^(^y(Tr,-n),ei}9(Ti) + ^(^'^(0), ei)6»(r„)^ , 

for each n G N“'", where t„ = nS constitutes a uniform grid of step size 5 > 0. Approximate 
solutions 9 can be constructed by solving for continuous piecewise linear (i.e. linear on 
each interval [T„,r„+i]) functions 9 satisfying 

= Xj,(0)(t„) = g(r„), forallneN+. 


(3.10) 
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As ly is a second order approximation of ly , it is not surprising that 0 is a second order 
approximation of 6. 

Lemma 3.8. Let {x,y) G 'KxY and g: M_|_ —> K piecewise with continuous second 
derivatives. If g{0) = 0, then there is a unique piecewise linear function 9 G (7(1^+) 


satisfying (3.10). Moreover, 6 is a second order approximation of the exact solution 9 of 
the Volterra equation Hy{9) = g in the sense that for each T G M+, 

sup \9{t) — 9{t)\ < C(5^, 

tG[0,T] 

where C is a constant depending only on T and g. 

The lemma will be used to show the numerical invertibility of the calibration operator 
'Hy{-,x). In this context, the smoothness assumption on the right-hand side g of the 
Volterra equation is satisfied if the yields are interpolated sufficiently smoothly. 

Proof. This follows from Theorem 3] and Example 2.4.5 and Theorems 2.4.5, 2.4.8], 
noting that the integral kernel ( |3.9[ ) is and strictly bounded away from zero along the 
diagonal by our assumption ('ll'(0),ei) = (A,ei) ^0. □ 


Note that solving for 9 can be performed efficiently because (3.10) is a linear system 
for {9{to), ..., 9{Tn)) of lower triangular form. 

3.9. Estimation of the affine coefficients. We first discuss how, in principle, the affine 
coefficients can be identified from covariations of yields and then present some practical 
considerations on the construction of estimators. 

Let r{t, t) denote the yield of a zero-coupon bond held from t to t + t, i.e. 


(3.11) 


r{t,r) = --\ogP{t,t + T), 

T 


for all t,T G K+. 


Then by Equation (3.5) the quadratic covariation of yields with maturity t, and Tj satisfies 


(3.12) 




(t) = - 'I’yipy 

TiTj 


E 






Assume that the left-hand side of (3.12) is known as a function of and Tj, and that 


the components of ipj, are functionally independent. Then the components of iky and the 
matrix Ay{X(t)) can be identified. The function diy usually determines the coefficients 
ay and j3y uniquely (see [Equation (3.4b) ). Furthermore, taking account of the admissi¬ 
bility conditions on the matrices Oy and a^ (see 15 Definition 2.6]) one can identify the 
IR+-valued components of X{t) and the matrix Oy. 


Note that (3.12) is derived solely from the diffusion coefficient of the yield dynamics 
and therefore is invariant under Girsanov’s change of measure. Thus, the coefficients 
ay,ay,/dy can be estimated from real world observations without specifying the market 
price of risk. Of course the market price of risk enters as a bias in the estimation, but 
the estimators do not depend on it. Moreover, under the model hypothesis the estimates 
do not depend on the choice of r^, Tj, which provides a means to reject ill-suited models. 

The remaining K-valued components of X(t) and the coefficient by do not appear in 
the quadratic covariations (3.12). We now discuss how they can be estimated, 
note that for one-factor models h 


First, 

is redundant and can be normalised to zero thanks 


to the Hull-White extension. In the multi-factor case only the first component of the 












12 PHILIPP HARMS, DAVID STEFANOVITS, JOSEF TEICHMANN, AND MARIO V. WUTHRICH 


vector by is redundant. Second, note that the short end of the forward rate curve gives 
a scalar condition on X{t), which allows one to fully identify X{t) if X{t) has only a 
single M-valued component. In the general multi-factor case, however, some components 
of by and X(t) remain undetermined. They may be calibrated to the prevailing market 
yield curve by regression methods. Alternatively, they may be estimated by econometric 
methods. However, these require a market price of risk specification. We do not discuss 


this topic here and refer to 22 for further details. 


In practise, the quadratic covariation (3.121 must be estimated from yields r{tmTi) 


given by the market for times = n5 and times to maturity t^. A naive estimator is the 
realised covariation, which is defined as 

n 

[r{-,n),r{-,Tj)]{tn) = ^ {r{tk,Ti) - r{tk-i,Ti)) {r{tk,Tj) - r{tk-i,Tj)). 
k=l 

Here also other estimators such as, e.g., Fourier estimators, as introduced by Paul Malliavin 


and Maria-Elvira Mancino, could be used, see, e.g., 12 for some recent developments. 


Fixing a time window of length M and a time one has 


(3.13) 


^n—M 

5 


— 'Sy{n)^ay'i>y{Tj) 


d 




k—1 m—n—M-\-l 

Therefore, for any time tn and any selection of times to maturity e stimators a 


f3y, X^{tn), ..., X‘^^{tn) can be constructed by solving for the best fit in Equation (3.13) 




4. Consistent recalibration of affine short rate models 

4.1. Overview. The constant parameter process y of the previous section is now replaced 
by a stochastic process Y = (F(t))t>o. The situation is particularly simple when Y is 
piecewise constant. In this case, the Hull-White extension is recalibrated to the prevailing 
yield curve (i.e., the yield curve given by the model with old parameters) each time the 
parameter process changes. Later on, the concepts are generalised to arbitrary parameter 
processes Y, resulting in our definition of general CRC models. Geometrically, these mod¬ 
els “locally look like” Hull-White extended affine short rate models with fixed parameter 
y. This is made precise in Section 4.8 A semigroup point of view is taken in [Section 4.10 


leading to an interpretation of CRC models with piecewise constant Y as splitting schemes 
with respect to a time grid for more general CRC models. 


4.2. Setup and notation. We recall the notion of admissible parameters from |Assump-| 
Eor all {s,x) G M+ x X and all admissible parameters (y,0) G V x C(U^), we 


tion 3.1 


let X = Xyg denote the unique solution on [s,oo) of the SDE (3.1) with 9{t) replaced by 


9(t — s) and initial condition X{s) = x. We assume that El is a Hilbert space of forward 
rate curves satisfying [Assumption 3.4| simultaneously for all y G Y. We fix a strictly 
increasing sequence of non-negative deterministic times {tn)neNo- 


4.3. CRC models with piecewise constant parameter process. We assume that the 
parameter process Y is piecewise constant, i.e. for each t G [f„, f„+i) we have Y{t) =Y{tn). 
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Definition 4.1 (CRC models with piecewise constant Y). A stochastic process {h,X,Y) 
with values in El x X x Y is called a CRC model if there exists a stochastic process 9 with 
values in (7(11+) such that the following conditions are satisfied, for each n € Nq: 

(i) The Hull-White extension 9 on [<„,<„+i] is determined by calibration to h{tn)'. 


h{tn){0) =i+ {X,X{tn)), 9{tn) = CY(t„) , X(t„)) , 


and for t G [<„,t„+i] 

9(t) = S{t - tn)9{tn). 

(ii) The evolution of X on t„+i] corresponds to the Hull-White extended affine model 
determined by the parameters (Y{tn),9{tn))'- 






t € [tyi , , 


where X^’n is the solution operator of SDE (3.1) defined in 


Assumption 3.1 is assumed to hold for the parameters {Y(t. 
(hi) The evolution of h on [t 


ij ^n+l 

Hull-White extended affine model: 


Section 4.2 Here, 
,0(t„)) G Y X (7(M+). 
is determined by X according to the prevailing 


h{t) —(0(t), Ar(t)) , t G [tn , tn+l] ■ 

We use the same symbols h and X as in [Equation (3.8) to denote CRC models. The 
abuse of notation is motivated by the fact that {h, X) evolves on , t„+i] according to (3.81 


with parameters (Y{tn),9{tn)). Note that the process X in Definition 4.1 is continuous 
because closed intervals [t„, t„+i] are used in point (ii). We emphasise that the recalibration 
step (|^ happens on a discrete time scale because the parameter process Y is constant on 
each [t„,t„+i]. By construction the process {h,X) is continuous at each time tn- 


4.4. Simulation. If we assume that a stochastic model for the evolution of the parameter 
process Y is specified, one can sample Y on the time grid (tn)neNo- Then CRC models 
as in Definition 4.1 can be simulated by applying iteratively steps (1^ 


Algorithm 4.2 (Simulation). Given h{to) and the process Y, calculate {h,X,Y,9) on 
the time grid (tn)nGNo by iteratively executing steps (|^-([^ of Definition 4.1 Abort with 
an error if the assumption in step (|ii| is not satisfied, for any n G Nq. 


The algorithm is illustrated in [Figure 4T] Note that the forward rate increments 
are calculated from increments of the affine factor process X, which can typically be 
simulated with high orders of accuracy and proper treatment of boundary conditions. 
These advantages are thanks to the affine structure of the CRC increments and are not 
available for general HJM models. 


4.5. Efficient updating of forward rate curves. Updating the curve of forward rates 
as prescribed by Definition 4. 11 involves calculating integrals on time intervals [0,t], for 
large values of r (see Section 3.7 for the formulas). A significant speed-up can be obtained 
when this update is done using the alternative formula provided by the following lemma, 
which involves only integrals over time intervals of length 6 = t„+i — tn- 
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{h{tn),X{tn+l),Y (i„), e{tn)) 




update h 


O [h{tn+l),X{tn+l),Y{tn),9{tn)) 


update Y 


{h{tn+l) , X (tn+l) ,Y (tn+l), d{tn)) 

Figure 4.1. Simulation of CRC models. Updating 9, X, h is done using 


(M, (1^ of Definition 4.1 respectively. Updating Y is done using 


the exogenously given model for Y. 


Lemma 4.3 (Efficient updating of forward rate curves). Definition 4-1 iii I can be rewritten 
as 


(4.1) h(tn+l) — S{S)h{tn) + 5((5)<i)y((^) - 

~ + J (^S{S — ds, 

where 5 = — t„. 

Proof. By conditions 0 and ( pn| of [Definition 4.1[ 

hit„+l)-S{S)hitn) = nYit„){S{6)9{tr,),X{tn+l)) - SiS)nY{t„){9{t„),Xitn)) 

= i - XY(t^) {S{5)9{t^)) - - {^'y^^,^yX{t^Yi)) 

— ^ + S{6)lY(t„) {9{tn)) + S{5)^'Yl^t„){S{5)'l!'Y(t^), X{tn)) ■ 

Now the assertion of the lemma follows from the relationships 

S{5)Iy{9)-Iy{S{5)9)= [ 9{s){S{S-s)-^'y,ei)ds, for all (5,0) G M+x C'(M+), 

Jo 

which can easily be verified from the definition. □ 


4.6. Bond prices and forward rates. 

Theorem 4.4 (Zero-coupon bond price and forward rate). Let {h,X,Y) be a CRC model 
as in \Definition i] with corresponding process 9. Define 

P{t, T) = e- h(t,s-t)ds^ ^ 0) = f -f (A, X{t)), B{t) = eJo 

Then the discounted price process 1 1 —>■ P{t, T)/B(t) is a ¥-local martingale, for each T > 0. 
In this sense, the bond market is free of arbitrage. Moreover, the following affine bond 
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pricing formulas hold: 

\og{P{t,T)) = -e{T - t) + [ e{t){s){'^^Yit){T-t-s),ei)ds 

Jo 

+ -t) + - t),X{t)), 

h{t,T) = i- [ 6l(f)(s)(^'y(4)(T-s),ei)ds-$y(()(r)-(«'V(t)(T),X(f)). 

Jo 

Note: the following proof shows the stronger statement that discounted bond prices 
are true martingales. 

Proof. On each interval [tmtn+i], the evolution of forward rate curves h{t) stems from 
a Hull-White extended affine short rate model. Therefore, for each T > 0, the discounted 
price process t >-)• P{t,T)/B{t) is a martingale on each interval Moreover, the 

process is continuously concatenated at the boundaries of the intervals. It follows that 
the process is a martingale on [0, oo). The affine bond pricing formulas are equivalent to 
h{t) = 'Hv(i)(0(t),N(t)), which holds by Definition 4.1 hi). 


□ 


4.7. Heath-Jarrow-Morton equation. 


Theorem 4.5 (HJM equation). Let {h,X,Y) be a GRC model as in Definition 4-1 with 
corresponding process 0 and assume that hit) € D{A), for each t > 0. Then the following 
properties hold: 

(i) the expression Cy{ t)(h(t )t Xit)) is well-defined and eguals Oft), for all t > 0; 

(ii) the parameters {Yit),9{t)) are admissible, for all t > 0; and 

(Hi) the process {h, X) is a strong solution of the following SPDE on El x X; 


(4.2) 


dh{t) = (^Ah{t)+f4l)^{X{t))yt + a^l^{X{t))dW{t), 


dXif) = y AY(t){Xit))dW(t) + {CY(t){h{t),X{t)){Q)ei+BY(t){X{t))yt. 
Proof. This follows from [Definition 4.1| and Theorem 3.5| 


□ 


4.8. Geometric interpretation. The consistent recalibration scheme has a nice geomet¬ 
ric interpretation. Forward rate curves of a Hull-White extended affine short rate model 
remain within the finite dimensional manifold with boundary given by 


f Oft + .s){^'yiT - s),ei)ds + i - ^yir) - (4'y(r),x) 
Jo 


(f, x) € M+ X 


as can be seen from [Theorem 3.^ These submanifolds foliate the space of forward rate 
curves or large portions thereof. Let h be a forward rate curve. Then, for every choice 
of functional characteristics {Fy,Ry), there is at most one leaf through h. However, if 
{Fy, Ry ) is allowed to vary, there are in general many leaves through h. A choice of leaf 
corresponds to a choice of foliation and thus to a choice functional characteristics {Fy, Ry). 

A CRC model is constructed by concatenating forward rate evolutions on leaves be¬ 
longing to different foliations. This allows the otherwise constant coefficients {Fy,Ry) to 
change over time. The result is an HJM model which is “tangent” to Hull-White extended 


affine short rate models. This is illustrated in Figure 4.2 
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Figure 4.2. Each affine short rate model foliates the space of forward rate 
curves into invariant leaves. CRC models are concatenations of forward 
rate evolutions belonging to different foliations or limits of such 
concatenations. 


4.9. CRC models. We generalise CRC models of [Section 4.3| to arbitrary parameter 
processes Y. In this section Y is not required to be piecewise constant as in the last 
sections. To characterise such models, we use the SPDE derived in [Theorem 4.5 


Definition 4.6 (CRC models). A CBC model is a process {h,X,Y) with values in HxXxY 
satisfying conditions of Theorem 4.5 with d{t) = Cy(t) {h{t), X{t)) for each t > 0. 


We may think of CRC models as continuous-time limits of the concatenations described 
Section 4.8| Note that these models satisfy all conclusions of [Theorem 4.4[ they 


are free of arbitrage because discounted bond prices P{t,T)/B(t) are local martingales 
thanks to the HJM drift condition 17 Theorem 6.1], and the short rate can be written 
as r{t) = h{t){0) = i + (A,X(t)). 


4.10. Semigroup interpretation. Assume that the parameter process Y is Markov on 



'P{t)f{ho,Xo,yo) = E [f{h{t),X{t),Y{t))\{h{0),X{0),Y{0)) = {ho,Xo,yo)] ■ 


Moreover, let Q = (Q(t))t>o denote the semigroup on x X x Y) obtained by holding 
the parameter process Y{t) = y fixed, i.e.. 


Q{t)f{ho,xo,yo) =E [f{h{t),X{t),yo)\{h{0),X{0)) = {ho,xo)] , 


where {h,X) are as in Theorem 3.5[ with y = y^. Finally, let TZ 
semigroup on x X x Y) describing the evolution of Y, i.e.. 


{TZ{t))t>o denote the 


Tl{t)f{ho, xo, yo) = E [/(/iq, xq, Y (t)) | Y(0) = yo] . 


Then, the concatenation {'R.{S)Q{S))^f of semigroups describes CRC models with a 
piecewise constant parameter process as in [Definition 4.1] Indeed, 


for all n S No, 


{TZ{S)Q{Sjrf{ho,xo,yo) = E [h(<„),X(t„),y(t„)], 






















CONSISTENT RECALIBRATION OF YIELD CURVE MODELS 


17 


where S = tn+i —t„ is the step size, {h{tn), X{tn)) is obtained by executing the simulation 
scheme of Section 4.4 and Y{tn) by sampling the Markov process Y on that time grid. 


4.11. Simulation of CRC models by splitting schemes. The semigroup interpreta¬ 
tion of [Section 4. 10| allows one to view Algorithm 4.2| as an exponential Euler splitting 
scheme for general CRC models as in [Definition 4.6 To see this, let/:IHIxXxY— 
be twice differentiable with derivatives being uniformly continuous on bounded sets and 
assume that El C dom(M). Then, under appropriate assumptions on Y, ltd’s formula holds 


for f{h{t),X{t),Y{t)) by 
of the generators 


131 

,G 


Theorem 4.17]. It follows that / lies in the common domain 
of the semigroups V, Q, TZ, and, if Y is independent of W, 

G'^f = G^f + G'^f. 


The exponential Euler splitting scheme with respect to this splitting is defined as 

VinS)/ « (exp((5^’'^)exp((5^®))"/ = (77.(5)Q((5))"/, for all n G Nq. 

By the considerations in [Section 4.10[ it coincides with the simulation scheme of |Section 4.4| 
The advantages of this simulation scheme in comparison to other methods are discussed 
in Sections (2.21 and [4.131 


4.12. Calibration of CRC models. In order to calibrate CRC models we need to esti¬ 
mate a time series of the parameter process Y from market data, and fit a model for this time 


series. Estimating a time series of the parameter process Y can be done as explained in Sec 


tion 3.9 for Hull-White extended affine models. The resulting time series Y consists of model 


parameters under a risk neutral measure, but they can be estimated from real world observa¬ 
tions since the estimators are obtained solely from the volatility of the forward rate process. 

Calibrating CRC models requires the additional task of selecting and fitting a model 
for the estimated time series of Y. This completes the model specification under a risk 
neutral probability measure. We do not discuss the market price of risk specification in 
this paper but refer to 


22 for more details. 


4.13. Robust calibration, consistency, and analytic tractability. It is time to ad¬ 
dress the question to what extent CRC models satisfy the interest rate modelling principles 
set forth in the introduction. First, we formalise the notion of consistency and the consistent 
recalibration property. We do this for CRC models, but it is unproblematic to generalise the 
definition to other forward rate models. Let prjj denote the projection of El x X x Y onto El. 

Definition 4.7. Let {h{t), X{t),Y{t))t>o be a CRC model and X C El x X x Y. Then 
the model is called consistent with I if {h{t), X{t),Y{t)) G I holds with probability one 
for any t > 0 and initial condition (/iq, xq, yo) G Y. Moreover, the model satisfies the 
eonsistent recalibration property with respect to I if the support of the law of h{t) contains 
prjj(X) for any t > 0 and initial condition (/iq, xq, yo) € I- 

The consistent recalibration property is equivalent to any open subset of pr][j(X) being 
reached at any time t > 0 with positive probability. Observe that the consistent recalibra¬ 
tion property does not hold on any reasonably large set I for Hull-White extended affine 
factor models as in [Section 3[ Indeed, as explained in [Section 4.8[ for any given initial 
value the process h remains within a finite dimensional submanifold of El. CRC models, 
on the other hand, enjoy the following properties: 
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Robust calibration: the robust calibration principle is satisfied perfectly. Indeed, 
the method described in [Section 4. 12| allows us to use the entire present and past 
market data of yields to select a model. Whenever possible, the parameters are 
estimated from realised covariations of yields, which allows one to bypass the 
usual inverse problems in calibration. 

Moreover, requiring parameters to remain constant throughout the life time 
of the model is less restrictive in CRC models than in the underlying affine factor 
models. The reason is that the parameters of the underlying affine model are 
turned into state variables of the CRC model. 

Consistency: the canonical state space I of CRC models is the subset of El x X x Y 
determined by the admissibility condition on the underlying Hull-White extended 


affine factor model (see Assumption 3.1). Under sensible specifications of the 
affine factor model, I n (H x {x} x {?/}) is large enough to contain all realistic 
market curves, for each fixed (x, ?/) S X x Y. If the Hilbert space El is continuously 
embedded in C^(M+), then X is also large in the topological sense of having 
non-empty interior. 

In this setup consistency holds by construction because the state process 
(h,X,Y) of CRC models does not leave the set X. The consistent recalibration 
property can be verified using standard arguments; the support of {h{t),X{t), Y{t)) 
is the closure of the reachable set of an associated control problem [^, the reach¬ 
able set is stable under the flows of the driving vector fields and their Lie brackets, 
and generically speaking, as soon as there is noise in the parameter process U, 
these vector fields together with their brackets span dense subspaces of El. The 
exact conditions are worked out for the Vasicek case in [Section 5.71 


Analytic tractability: the simulation scheme for CRC models (Algorithm 4.2) 


transfers the task of sampling state variable increments to a finite-dimensional 
setting. Namely, instead of simulating forward rate increments from an infinite¬ 
dimensional space, it is sufficient to simulate increments of the finite-dimensional 
processes X and Y. This allows one to take advantage of the existing high-order 
schemes for the simulation of affine processes. Note that all operations in AU 
jgorithm 4.2| which involve infinite-dimensional objects are deterministic. The 
complexity for simulation reduces dramatically when high order methods for 
finite dimensional equations are applied which are often not at hand for infinite 
dimensional equations. Additionally often exact schemes are available in the affine 
finite dimensional setting, e.g., for CIR or Wishart type processes besides of course 
Gaussian processes. 


5. Consistent recalibration of Vasicek models 

5.1. Overview. We describe CRC models based on the Hull-White extended Vasicek 
model in full detail. Moreover, we show using semigroup theory that the simulation scheme 
of Section 4.4 converges to the CRC model of Definition 4.6|in the continuous-time limit. 


5.2. Setup and notation. We use the setup of Sections 3^ and |4.2[ setting X = K, 
£ = 0, A = 1. We do not specify the parameter space Y, yet, but we assume that for 
each {x,y) e X x Y, the volatility and drift coefficients are given by Ay{x) = Uy G [0,oo) 
and By{x) = PyX with Py G (—oo, 0). For simplicity, we choose equidistant grids of times 
tn = n6 and times to maturity r„ = nS, for all n G No, where <5 is a positive constant. 
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5.3. Hull-White extended Vasicek models. For each parameter {y,6) S Y x C(M+), 
the SDE for the short rate process is 


(5.1) 


dr(t) = {0{t) -I- I3yr(t))dt + y^dWit), 


where W is one-dimensional (J^(t))t>o-Brownian motion. Assumption 3.1 is satisfied for 
each parameter {y,0). The functional characteristics {F, R) from Section 3.4 are 


^ Ryi.u) = Pyu, 
and the solutions (ihyjdiy) of the corresponding Riccati equations are 
<^>y{t) = ^ {2Pyt - + 3 + ^ (l - , for alH > 0. 

Py 


for all M G 


By Theorem 3.3 the forward rates in the Hull-White extended Vasicek model (5.1) with 
fixed parameters {y,0) are given by h{t) = 'Hy{S{t)9,r{t)) G (7^(11+), where 

nyie, x){t) = r ^ (l - 

Jo ‘^Py 

for all {x,9,t) G M x C'(K+) x IR_|_. Due to the simple structure of the integral kernel 
there is a closed-form expression for the calibration operator, 

(5.2) Cy{h){T) = h'{T)-Pyh{T)-:^{l-e‘^^y^), for all (h, r) G C^(M+) x K+. 

2Py 

This can be verified using the definitions. Note that the calibration operator does not 
depend on x. Therefore, we dropped x from the notation Cy{h,x). 

The HJM drift and volatility from [Section 3.6 


are 


(5.3) /i™(r) = -^e^“"(l-e'^“"), a™(T) = for all r G K+. 

Py 

Note that these expressions do not depend on x, which is why we again dropped x from 
the previous notation •^'^(x)(r),CT^'^^(x)(r). The HJM equation for forward rates then 
reads as 


(5.4) 


dh{t) = ^Ah{t) + p. 


HJM 

y 


dt 


HJM 

y 


dW{t). 


5.4. Vasicek CRC models. Since the factor process is a function of the forward rate 
process, i.e., X{t) = r{t) = h{t,0), the corresponding CRC models can be characterised 
by the process {h,Y) instead of {h,X,Y). Thus, in accordance with Theorem 4.5 and 


Definition 4.6 a process (h, Y) with values in IH x Y may be called a CRC model if h 


satisfies the SPDE 

(5.5) dh(t) = (^Ah{t) + + CTp™dlE(t), 


with drift Pylt) volatility up™ defined in (5.3). Beyond the obvious requirement that 
these quantities are well-defined, for all t G K+, no further conditions are needed. In other 
words, the maximally admissible set I in the Vasicek case is the entire Hilbert space H. 
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5.5. Simulation of Vasicek CRC models. Given the parameter process Y with values 
in Y, the CRC model is simulated as described in [Algorithm 4.2| The following obser¬ 
vations make the algorithm particularly effective. First, the state process A is a function 
of the forward rate and can be eliminated as a state variable. Second, the short rate 
process can be simulated exactly. Indeed, in the model with constant parameter y, r(t) 
is normally distributed, 

r(t) Affe^^^ro -h [ e^y^^~^'^9{s)ds,- l) 

\ Jo ^Py 


Third, inverting the Volterra integral operator can be avoided by using closed-form 
expression (5.2) of the calibration operator. 

Discretisation is done on the uniform grid = r„ = 6n for a choice of finitely many 
times to maturity r„. Integrals are approximated to second order by the trapezoid rule, 
which leads to a global error of order one (see Section 5.6 and Section 7.7). The resulting 
scheme works as follows. 


Algorithm 5.1 (Simulation). Given (/i(0), Ah(0)) and the parameter process Y, execute 
iteratively the following steps, for each u G Nq: 

(i) The values of 0(t„) = CY(t„){h{tn)) at times to maturity 0 and 5 are calculated using 

0itn)iO) = Ah{tn)iO) - l3Y{t^)Htn){0), 

= Ah{t^){5) - pYit^)h{t^m - (1 - , 

and lY(t„)(Q{tn)){5) is approximated by the trapezoid rule as follows: 

^y(t„)(^(in))(<5) = (e^"<‘">'0(tn)(O) + 0(t„)(5)) . 


(ii) A sample r(t„+i) is drawn such that conditionally on A(f„), r{tn+i) has normal 
distribution 



(iii) (h(t„_|_i), A/i(t„+i)) is calculated from (/i(t„), Ah(<„),r(t„+i)) using 


Lemma 4.3 


h(t„+i)(r) =/i(t„)(5-f r)-f - (l - ^ 

+ -f r(t„+i) +lY(t„){0{tn)){S)^ , 

Ah{tn+l){T) = Ah(ti){S + t) 

_l_ ^Y{tn) 7g/3-K(t^)T _|_ g2/3y(t^)(T-|-5) _ g2/3y(t^)T _ g/5i'(t,i)(<5-|-T)'\ 

f^Y{t„) ^ 

+-f r(t„+i) -f ( 6 »(t„)) (5)^ . 

Here, h{tn+i) must be calculated at all times to maturity t^, whereas Ah{tn+i) is 
needed only at tq = 0 and ti = <5. 
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5.6. Convergence of the simulation scheme. In this section, we show that scheme of 


We are not aiming for the highest generality. Instead, we show how the results follow from 
standard semigroup theory. 


Section 5.5 converges to the CRC model (5.5) as the size <5 of the time grid tends to zero. 


Assumption 5.2. The parameter process Y takes values in Y = and satisfies 
(5.6) dV(t) = (AV(t) + ^(Y(t)))dt + a(Y(t))dW(t), 

where A : —>■ is a linear mapping generating a semigroup of contractions on fx € 
Cl^(R^;RP), a € and W is q-dimensional J'(t)-Brownian motion, indepen¬ 

dent of W. We write for bounded functions with bounded derivatives of all orders. The 
above SDE has a unique solution for any initial condition Y(s) = y, where (s, y) G K+ x Y. 


Assumption 5.3. The mappings y i—?> and y ^ Py are of class C^iW’) and 

supj^gY /3y < 0 holds. 


As Vasicek CRC models can be characterised in terms of {h, Y) instead of {h, X, Y), the 
semigroups V, Q, and TZ from Section 4.10 are now assumed to be defined on Ch(IHI x Y) 
instead of C;,)!! x X x Y). Recall that V describes the joint evolution of the process {h, Y), 
Q the evolution of h with Y fixed, and 71 the evolution of Y with h fixed. 


Theorem 5.4. There exists a separable Hilbert space El of continuous functions on K_|_ 
and a Banach space B of continuous functions on Hx Y such that 7^, Q, and TZ are strongly 
continuous semigroups onB. Moreover, for eachT G K_|_ there exists a constant C such that 

\\V{t)f - (7^(^/n)Q(^/n))”/llB < ll/IU' > V all fGM',tG [0,T],n G N+, 
where B' is a Banach space which is densely and continuously embedded in B. 

The space B' is large enough to be relevant in applications: any function on Ho x Y 
belongs to B', where Hq D H is defined in the proof below. 

Proof. We proceed as in 14 and 16 . Let be a strictly increasing sequence of real 

numbers strictly greater than 3. For each i € No, define a separable Hilbert space H^ by 

Hi = < /i e Lioc^ G ^loc a-nd [ h^^\T)‘^{l + < oo, Vj = 1,...,i + 1 > , 

[ J{0,oo) j 

where denotes the space of locally integrable functions on (0, oo). Every function in 
Ho is continuous, bounded and has a well-defined limit h(oo) = limT-_ioo h^r). The scalar 
product on Hi is defined by 

= ^i(oo)h2(oo) + ^ [ h[^\T)h':^^\T){l + ryMr. 

d(O.oo) 

For each C > 0 and fc, z S No, we define the space B^,(Hi x Y) as the closure of C^(Hi x Y) 
under the norm 

k 


Bi(H.xY)=Zl (cOs1i(CI|/j||hJ 


III/IIy) ll^■’7(^)^/)llL((HiXY)3)• 


Together with [Assumption 5.2 and 5.3 this implies that the conditions of 14 Sections 
3.1.1 and 3.1.2] are satisfied for SPDE (5.5), (5.6) characterising the evolution of {h,Y). 
(Note that fdy needs to be bounded away from zero for and to be bounded 
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with bounded derivatives.) Thus, this SPDE admits unique solutions on each space Hi x Y, 
given that the initial condition is smooth enough. The same applies to the SPDE for h 
with fixed y and the SDE for Y with fixed h. 

Eix (^0 > C > 0 and define H = H2, B = Bg“(H2 x Y), and B' = B4(Ho x Y). Then 
i'P{t))t>o, {Q{t))t>o, and {TZ{t))t>o are strongly continuous semigroups on B by 14 


14 


Lemma 13] and quasicontractive by |14[ Lemma 7]. Their generators are denoted by 
, Q^, and By the same lemma, B' is stable under {V{t))t>o- Together with 
Theorem 11] this implies that for each / € B', the expressions 

g^r{t)f, g^v{t)f, g^g^r{t)f, g^g^vit)/, g^g^vit)/ 

are well-defined with B-norm bounded uniformly in t € [0,r] and f = Q^f + g^f. 
Thus, the splitting is of formal order one and the result follows from [20[ Theorem 2.3 
and Section 4.4]. □ 

5.7. Consistent recalibration property. If the coefficient f3 in the HJM volatility 
is stochastic, one would expect the forward rate process to reach every point in the 
Hilbert space with positive probability, i.e., the consistent recalibration property holds. 
This is made precise here. 

Assumption 5.5. For each initial condition y(0) = j/o S Y and each T > 0, the support 
of Yt is all of Y. Moreover, {f3y : y gY} contains an interval [f3, 00) for some /3. 


Assumption 5.6. The Hilbert space 
finite abscissa 


is contained in LL(M+), and each /i S H has a 


abs(/i) = inf ^^(3 gR : J h{T)e^'^dT < oo| < 00. 

The condition in [Assumption STb is mild; it is satisfied by the weighted Sobolev spaces 


16 and in particular by the space H of Theorem 5.4 The above assumptions imply 


the consistent recalibration property, as the following theorem shows. 

Theorem 5.7. The consistent recalibration property is satisfied for the Vasicek CRC 
model (5.5), (|5.6|) with respect to the state space I = H x Y. 


Proof. Let {h,Y) be the solution of (5.5) and (5.61 with initial value (/iq, Tq) and let T > 0. 


By 27 Theorem 1.1] the support of {hx, Yt) is equal to the closure Lt of Ct, where Lt is 


the reachable set at time T of the following control problem: in (5.5) and (5.6) the Brownian 


motions are replaced by piecewise continuously differentiable control functions. Let (h,Y) 
be the solution of (5.5) and (5.6) for the zero control. Taking variations in the control for 


Y implies that {hx} x Y C Ct thanks to Assumption 5.5 Adding variations in the control 
of h improves this to ({/it}+ span{ : y G Y}) x Y C Ct- The set spanjcTy ■. y gY} 
is dense in H because its orthogonal complement vanishes by Proposition 1.7.2] and As¬ 


sumption 5.6 Therefore, H x Y C Ct, and the consistent recalibration property holds. □ 


5.8. Example. We present an example of a Vasicek CRC model based on 10 . In this 


model, the volatility is stochastic, but the speed of mean reversion is not. Therefore, 
the conditions of [Theorem 5.7| are not satisfied, and it turns out that the model admits 
a finite-dimensional realisation. The explicit formula for bond prices in this model will 
serve as a reference for showing convergence of the numerical simulation scheme for CRC 
models in the continuous-time limit. 
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The parameter process F is a CIR process with values in Y = K+ given by SDE 
( |5.6[ ) with a possibly correlated Brownian motion W and coefficients = m + and 
a{y) = ay for some m > 0, fj, < 0, and cr > 0. The Vasicek drift and volatility in the HJM 
equation are given by Equation (5.3) with /3y = /3 for some constant /3 < 0 and Uy = y, i.e., 


HJM 


(r) = -|e^"(l-e''"), a™(r) = for all r e M+. 


If h{0) S there is a closed-form solution of CRC equation (5.51, 

h{t,T) = h{0,t + T)- [ ds-k [ ^/Y{^el^'^^+*-‘UW{s). 

Jo P ^ ^ Jo 

Setting ^(t) = /p and r{t) = h{t,0), this can be rewritten as 

(5.7) h{t, t) = /i(0, t -I- t) -I- (r(t) - h(0, t)) + ^ - e^'^) ^(t). 

Setting T = 0 in [Equation (5.5)| and plugging in [Equation (5.7) yields 

dr(t) = (Ah(0,t) — j3h{Q,t) + + C(^)) + \/Y{t)dW{t). 

Summing up, the process X = (r, Y) is given by the SDE 

dr{t) = ^^(0,t) - I3h{0,t) + I3r{t) + ^(t)^ dt + y'Y{t)dW{t), 

dm = {m{t)+Y{t))dt, 

dY{t) = (m + yY{t))dt + ayjY{t)dW{t), 

where h{Q) G C^(U+) is a given initial forwar d rate curv e. It follows that X is an affine 
factor process for the short rate as described in Section 3 with d = 3, £ = 0, A = (1, 0, 0)^. 
Thus, the CRC model has a finite-dimensional realisation. If cr = 0, the affine bond pricing 
formula is particularly simple: bond prices are given by 

J'j _ g/t^(e'^'''“*’ft(0,t)-Ii(0.s))ds+/3“^(l-e'^*)r(i)-i/3“^(l-e'^*)^{(t)^ 

where ^{t) is deterministic and satisfies 


r(o) ~ I ~^~ if^<0 


(5.8) 


m = 


r(o) 


2/3-/X ' 2/3/x(2/3 -/x) 

g2/3i jji — 2j3t — l) 


2/3 4^2 

and r(t) is normally distributed with mean 


if /X = 0, 


(5.9) 

and variance 


e^*r{0) + f (^Ah{0, s) — Ph{0, s) + ^{s))ds, 

Jo 




(5.10) 


2/3 
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5.9. Calibration of Vasicek CRC models. As outlined in [Section 4.12[ we first con¬ 
sider y as fixed and suppress the dependence on y in the notation. For any selection of 


times to maturity Ti,Tj, estimators for a,/3 can be obtained as described in Section 3.9 


by solving for those a, P which achieve the best fit in (3.13), i.e. 

(5.11) 


^n. in. — 


M 


- 1 - 1 

Pn pTj 


Pn + P^Tf/2 PTj + P^Tf/2 


Pn 


PTj 


Varying the calibration time tn creates a time series of coefficients a(t„),/3(t„) for which 
we need to specify and calibrate a model. Some models are described in [Section 7.5[ below. 

6. Consistent recalibration of Cox-Ingersoll-Ross models 
6.1. Overview. We give a brief overview of CRC models based on CIR short rates. 


The overview is sufficient to set the notation for the empirics in Section 7 A detailed 


description is provided in the online appendix to this paper. For comparison, we briefly 
digress to the CIR-I—f model and its CRC version. 

6.2. Hull-White extended Cox-Ingersoll-Ross models. We use the setup of Sec¬ 
tions 3.2 and [4.2[ setting X = IR+, £ = 0, A = 1. We do not specify the parameter space 
Y, yet, but we assume that for each {x,y) S X x Y, the volatility and drift coefficients 
are given by Ay{x) = UyX and By{x) = PyX for some ay G (0,oo) and Py G (—oo,0). 
For simplicity, we again choose equidistant grids of times = n5 and times to maturity 
Tn = nS, for all n G Nq, where 5 is a positive constant. 

The CRC algorithm is similar to the Vasicek model, with the following important 
differences: 


The admissibility condition in Assumption 3.1 is satisfied if and only if 0(t) > 0, 
for all t G K+. This condition can be problematic in practise, as discussed in 


Section 7.4 Moreover, if this condition is expressed in terms of forward rates h(t) = 


'Hy{0{t), X{t)) instead of Hull-White extensions 0(t), it becomes apparent that 
the set of admissible forward rate curves depends on the parameter y. This makes 
it difficult to formulate convergence results similar to those in the Vasicek case. 

• In contrast to the Vasicek model, there does not seem to be a closed-form ex¬ 
pression for 0 = Cy{h,x) because the Volterra kernel in the integral operator Hy 
is more complicated. Therefore, the Volterra equation is solved by numerical 
approximation of order two using a discretisation in the time to maturity. 

6.3. CIR-I—h models in the CRC framework. In the CIR-I—I- model Section 3.9], 
also known as deterministic shift-extended CIR model, the short rate process is defined 
by r{t) = X{t) + 0{t), where A is a CIR process and 0 is a deterministic function of time. 
Note that this is a different time-inhomogeneity than the one described in [Section 6.2[ 
In particular, the factor process X is time-homogeneous and does not coincide with the 
short rate. The HJM equation of the CIR-I—I- model is 


dh{t) = [Ah{t) 


M™(A(t))) dt + a™{X{t))dW{t), 


dX{t) = {by -f PyX{t)) dt + \Ja^l^X{^dW{t), 


(6.1) 
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where and are the same as in the CIR case. In the CRC extension of this 

model, y is replaced by a stochastic process Y. This model has both advantages and 
disadvantages over the consistently recalibrated CIR model: 


• The SDE for X does not depend on h. Therefore, existence and uniqueness of X 
can be shown by standard methods. Then a mild solution h can be constructed 

Section 6.1]: 


by stochastic convolution 13 


h{t)=S{t)h{0)+ St-sti^l^{X{s))ds+ St-sa^l^{X{s))dWis). 


• The function 9 is allowed to assume negative values and can be calibrated to a 
given yield curve without having to invert a Volterra integral operator. However, 
this calibration requires one to know the current value of X. This can be seen 
from the equation for forward rate curves 


h{t) = S(t)9 - by'l’y - %X{t), 


where 'i'y is the same as in the CIR case (c.f. Section 3.4). In contrast to the CIR 
model, the process X is not directly observable from the short end of the term struc¬ 
ture. Therefore, X{t) and the parameters ay and fdy have to be estimated jointly 
from realised covariations of yields as described in |Section 3.7[ Moreover, in con¬ 
trast to the CIR model the parameter by is not redundant and has to be estimated 
using the same methodology as for general multi-factor models (see Section 3.7). 


7. Empirical results 

7.1. Overview. CRC models based on Vasicek and CIR short rates are calibrated to Euro 
area yield curves. Properties of the calibrated models are studied in comparison to market 
data and models without consistent recalibration. Our empirical findings show that the 
assumption of constant parameters in the Vasicek and CIR models is too restrictive. There¬ 
fore, the additional flexibility provided by CRC models is a useful tool to better capture 
the market dynamics. This is also reflected in better fits of the covariance matrix of yields. 


7.2. Description of the data. We consider the zero-coupon yield curves released by the 
European Central Bank (ECB) on a daily basis. The yields are estimated from AAA-rated 
(Fitch Ratings) Euro area central government bonds being actively traded on the market. 
Estimation is done by the ECB using the Svensson family of curves, see 31 33 . Data 
is available from September 6, 2004, and we set April 1, 2014 to be the last observation 
date. In total, this results in 2454 observed yield curves with times to maturity ranging 


from 3 months up to 30 years. Yields are continuously compounded (c.f. Equation (3.11)) 


and denoted by r(t,r), with r being the time to maturity. A selection of yields is shown 
in Figures [TT] and [7?^ The short rate is approximated by the yield with the lowest time 


to maturity (3 months) and is depicted in Figure 7.1 


7.3. Calibration of CRC models. The CRC models based on Vasicek and CIR short 
rates are calibrated as described in [Section 5.9| and the online appendix. Time steps 
S = 240“^ of one business day and time windows of M = 100 business days are used for the 
calibration. The choice M = 100 is a compromise between accuracy and over-smoothing 


and gives reasonable results over the time horizon of roughly a decade (see Figures 7.3 and 
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Zero-coupon bond yields (%) 



Figure 7.1. Historical zero-coupon yields estimated by the ECB for 
various times to maturity from 06/09/2004 to 01/04/2014. We use the 
3-month yields (r = 0.25) as a proxy for the short rate. 


Zero-coupon bond yields (%) 



Figure 7.2. Zero-coupon yield curves estimated by the ECB for various 
observation dates. 


7.4). For Ti ^ 1 one immediately obtains from Equation (5.11) and its CIR counterpart 


the estimator 

(7.1) a{t) = 

in the Vasicek case, and 


[^(LTi),r(-,ri)](f) - [r{-,Ti),r{-,Ti)]{t- SM) 
■ "" 6M 


a{t) 


[r(-,ri),f(-,Ti)](t) - [f(-,Ti),f(-,ri)](t - 8M) 
^ Em=0 h) 


(7.2) 
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Volatility in the Vasicek model (%) 



time (t) 


Figure 7.3. Volatility parameter of the Vasicek model estimated by 
using a time window of M = 100 yields with time to maturity ti. 


in the CIR case, where the quadratic variation is estimated by (5.111. On the other hand 


taking r 2 ^ 1, one can solve (5.111 and its CIR counterpart for /3 and obtain the estimator 
(7 3 ) Rif) = _i_ _ SMa{t) _ 

in the Vasicek case, and 


/3(i) = ( [H-,T-2) ,H-,T2)]{t) - [r{-,T2),r{-,T2)]{t- SM) 


(7.4) 


ST,^=o - Sm,Ti) 


\/W) ( [^(LT-2),T(-,r2)](t) - [r(-,r2),f(-,T2)](t- 5M) 


T2 


Em=o '^i) 


in the CIR case. The resulting trajectories of the estimated Vasicek volatility 
CIR volatility 


/a are shown in Figures 7.3 and 7.4 


o and 

The trajectories of the speed of mean 


reversion —/3 for both models are plotted in Figures [71^ and [7161 

It turns out that for most parts of the data, a and a do not depend too much on the 
times to maturity used in the estimation, whereas j3 does. Typically, smaller times to 
maturity result in larger values of —/3, as shown in Figures [7)5] and 7.6 This means that 


one-factor Vasicek and CIR models are not flexible enough to reproduce the ECB yield 
curve movements in full accuracy, and so the choice of times to maturity used in the 
calibration procedure may have an impact on the results. In particular we set ti = 0.25 
and T 2 = 2 (i.e. 3 months and 2 years). 

The dependence of the estimated parameter /3 on the choice of times to maturity 
suggests to use multi-factor models as building blocks for CRC models. In the empirical 
part of this paper, we aim for a detailed understanding of the one-factor case and leave 
the extension to multiple factors for future research. As pointed out, accurate modelling 
of long-term rates might require one to enlarge the model class even further to allow the 
forward rate volatilities to decay slower than exponentially. 
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Volatility in the CIR model (%) 



Figure 7.4. Volatility parameter of the CIR model estimated by (7.2) 
using a time window oi M = 100 yields with time to maturity ti. 


Speed of mean-reversion in the Vasicek model 



Figure 7.5. Speed of mean-reversion parameter of the Vasicek model 
estimated by (7.3) using a time window of M = 100 yields with times 
to maturity ti = 0.25 and various values of T 2 . 


By the theory of Hull-White extensions, an exact match to the initial yield curve is 
always achieved, but the corresponding time-homogeneous models often do not match the 
initial yield curve well (Figures 7.7 and 7.8). This is not surprising as they are calibrated 
to the yield curve dynamics and not to the initial yield curve. This separation of dynamics 
and initial calibration is actually one of the strengths of our approach. 

To model the dynamics of the Vasicek coefficients a, (3 and the CIR coefficients a, /3, we 
use geometric Brownian motions and/or CIR processes, as laid out in Section 7.5 Note that 
the assumption of constant coefficients, which is implicit in affine factor models without 


the CRC extension, is not realistic over long time horizons in view of Figures 7.3 and 7.4 
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Speed of mean-reversion in the CIR model 



Figure 7.6. Speed of mean-reversion parameter of the CIR model 
estimated by (7.4) using a time window of M = 100 yields with times 
to maturity ti = 0.25 and various values of T 2 - 


Zero-coupon yields (%) 



time to maturity (x) 


Figure 7.7. Calibrations of some homogeneous and Hull-White extended 
models as of 1 April 2014. Vasicek 1 and CIR 1 are homogeneous models 
calibrated to the yield curve dynamics using (7.1|-(7.4) with n = 0.25 
and T 2 = 2. Vasicek 2 and CIR 2 are homogeneous models calibrated 
to the prevailing yield curve by least squares. The Hull-White extended 
models match the initial yield curve exactly. 


7.4. Negative levels of mean reversion. A problematic aspect is that the time- 
dependent drift 9 obtained by the calibration to the initial yield curve can assume 
negative values, which are not admissible in the CIR model. The problem occurs mostly in 
low interest rate scenarios with partially inverted yield curves (Figures 7.7 7.8 and 7.9). 
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Time-dependent drift (%) 



Figure 7.8. The Hull-White extensions 9{t) corresponding to the 
models of Figure 7.7 Note that 9{t) assumes negative values, which is 
typical in situations where the yield curve is inverted at the short end. 


initiai time-dependent drift in the CIR modei (%) 



Figure 7.9. The historical values of 9{t){0) in the CIR model calculated 
using the estimates of (3 in Figure [7?6l Negative values occur frequently 
in 2009 and 2012 to 2014. They are problematic for reasons laid out in 
ISection 7Hl 


In contrast, negative values of 9 are allowed in the Vasicek model and might even be 
desirable for modelling bond markets with negative interest rates. 

As only the short (left) end of 9 is relevant for CRC models, at each step of the 
simulation scheme (cf. Algorithm 5.1), it is sufficient to understand how 0(0) depends 
on the prevailing forward rate curve and the coefficients of the affine factor process. The 
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general formula for 0(0) is 

e{Q)=Cy{h,x){Q) = 


ei' 


{h\0)-F'{0)-X-{R'y{0)-X,x)) 


which follows by differentiating the relation h = 'Hy{9,x) with respect to the time to 
maturity r and evaluating at r = 0. In the Vasicek and CIR case, the formula becomes 
0(0) = h'{0) — Pyh{0). This shows that the problem can be alleviated to some extent by 
artificially choosing higher levels of —/3, resulting in higher values of 0(0). For this reason 
we set T 2 = 2 instead of higher values of T 2 in the calibration of (3 (Figures 7.5 and 7.6). 
Despite this correction, 0(0) remains negative for most parts of 2009 and 2012-2014. 
During these periods of time, the CIR model cannot be fitted simultaneously to the yield 
curves as well as their volatilities and has to be rejected. 


7.5. Models for parameter evolutions. There are very few restrictions on the choice 
of parameter process. It can be chosen exogenously for scenario based simulation or 
calibrated to the market, and it is not restricted by the fit to the initial term structure. 

We consider here four models for the evolution of the Vasicek and CIR parameters: a 
reference model where the parameters are constant, two toy models with constant mean 
reversion and time-varying volatility, and one fully stochastic model which is calibrated 
to the market. In the Vasicek case, the four models are: 

(VI) a Hull-White extended Vasicek model with constant coefficients jSy = /Iq and ay = oq; 
(V2) a Vasicek CRC model with constant mean reversion coefficient Py = /3o and deter¬ 
ministically increasing volatility given by ay = aoj/, V(t) = 1 + 3t; 

(V3) a Vasicek CRC model with constant mean reversion coefficient f3y = /3o and stochasti¬ 
cally increasing volatility given by Uy = y, dY{t) = (4ao — Y{t))dt + a^yY{t)dW(t), 
V(0) = tto, cr = 3 • 10“^; and 

(V4) a Vasicek CRC model with stochastic coefficients given by geometric Brownian 
motion: = yi, ay = 1 / 2 , dYi = yiYi{t) + aiYi{t)dWi{t), dY 2 = H 2 Y 2 {t) + 

cr 2 Y 2 {t)dW 2 {t). The coefficients ^ 1^2 and CTi ^2 are deterministic and calibrated to M 
observations as described in [Section 7.31 


Note that models (\^ (\|^ both describe scenarios where the standard deviation of 
the noise in the short rate process doubles over the period of a year. Indeed, in (\[^ the 
volatility coefficient a increases to four times its initial value, and in (V[^ the level of 
mean reversion of a increases to four times its initial value. 

Models (\§ and (V[^ are special cases of |Section 5.8[ In (\|^, which corresponds 
to TO = 3ao, y = 0, and tr = 0, there is an explicit formula for the moment generating 
function of the short rate process. By equations ( |5.9[ ) and ( |5.10| ), it is given by 

g/3o(‘-=) (^hg{s)-l3hois)+^{s))ds+^^{t) 


(7.5) E 




_ ^ef’o^rjro+vfo 


Tj € ]R, t € 




where 


^(t) = ao 


e2/3oi _ 1 Sag (e2/3ot _ 2j3ot - l) 


+ 


2/3q 4/?q 

Model (\|^ corresponds to Section 5.8 with to = 4ao, /i = — 1 and cr = 3 • 10“^. 

The semigroup approach of Theorem 5.4| implies convergence of the simulation scheme 
for (\§]) and for (\Q with Yi replaced by Yi -|- e for some e > 0. In our numerical sim¬ 
ulations, we observe convergence for all models (see Section 7.7), including the following 
CIR counterparts of the models just described: 
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(CIRl) a Hull-White extended Cox-Ingersoll-Ross model with constant coefficients Py = Pq 
and ay = ao\ 

(CIR2) a Cox-Ingersoll-Ross CRC model with constant mean reversion coefficient Py = Po 
and deterministically increasing volatility given by ay = agy, R(t) = 1 -I- 3t; 

(CIR3) a Cox-Ingersoll-Ross CRC model with constant mean reversion coefficient Py = Po 
and stochastically increasing volatility given by ay = y, dY{t) = (4ao — Y{t))dt-\- 
asjY{t)dW{t), R(0) = ao, cr = 5 • 10“^; and 

(CIR4) a Cox-Ingersoll-Ross CRC model with stochastic coefficients given by geomet¬ 
ric Brownian motion: Py = j/i, ay = y 2 , dYi = /iiYi(t) -|- cfiYi(t)dWi{t), 
dY 2 = y, 2 Y 2 {t) + cr 2 Y 2 {t)dW 2 {t). The coefficients /ri ,2 and cri ^2 are determin¬ 
istic and calibrated to M observations as described in ISection 7.31 


7.6. Implementation. Simulating CRC models requires iterative sampling of the under¬ 
lying affine short rate process and recalibrating Hull-White extensions. This is explained in 
detail in |Algorithm 5.1| for the Vasicek model and in the online appendix for the CIR model. 
The algorithms can be parallelised on a path-by-path level. Parallelisation on lower levels 
does not pay off because the individual time steps are dependent on each other. In our im¬ 
plementation in R, generating 10^ paths with 240 time steps on a cluster of 48 times 2.2GHz 
processors takes around 10 minutes in the Vasicek case and 20 minutes in the CIR case. 


7.7. Convergence analysis. Theorem 5.4 predicts first order convergence of the simu¬ 
lation scheme under suitable assumptions on the model. The objective of this section is to 
demonstrate this convergence in numerical examples for the models described in |Section 7.5| 
The simulation is started with the initial forward rate curve hg of September 2, 2013. 
The parameters /3o, og, og, 0o are calibrated as in [Section 7.3| with a time window of 
M = 100 observations. Then the moment generating function of the short rate r(l) 
after one year is calculated by Monte Carlo simulation with 10® sample paths. In the 
model (\|^, the exact value of the moment generating function is known and given by 


Equation (7.5) In the other models, a reference value is calculated by extrapolation from 


the Monte Carlo estimates. The resulting errors are shown in Figures [7. 10[ |7.li] and |7.12| 
As expected from Theorem 5.4| we observe first order convergence for models (\|^ and 
(\Q. The errors in Figure 7.12 indicate convergence also in the CIR counterparts. 

7.8. Distributional properties. Making parameters in the CIR and Vasicek model 
stochastic in the sense of CRC models has considerable impact on the distribution of short 
rates and prices. As an example, statistics of the short rate r(l) obtained by simulation 
are presented in Table 7.1 In the mod els (A[T]) and (\|^ with deterministic parameters, 
the short rate process is Gaussian (see Section 5)^ . As expected, the simulations show 
skewness and excess kurtosis values close to zero. In contrast, leptokurtosis appears in the 
models (\|^ and (\Q with stochastic parameters. In the CIR examples, the distribution 
of r(l) is also affected considerably by the stochasticity of the parameters. 


7.9. Covariation of yields. A further example where empirical differences between 
CRC and non-CRC models become apparent is the covariation matrix of yields. By the 
Dubins-Schwarz theorem, the covariation process determines the distribution of yields and 
therefore the prices of bond options. On short time intervals, the covariations are closely 
related to covariances, which are a key determinant of the prices of call and put options. 

In the Hull-White Vasicek and CIR models without the CRC extension, the covariation 
matrix of yields with different times to maturity has rank 1. This is in stark contrast to 
























log. abs. error log. abs. error 


CONSISTENT RECALIBRATION OF YIELD CURVE MODELS 


33 


Absolute error of the Monte Carlo MGF for model (V2) (log. scales) 



Figure 7.10. Absolute error (log-log plot) of the Monte Carlo estimate 
of the moment generating function E for model (’Vj^. This is 

calculated as the absolute difference between the estimate and (7.5) for 
different values of 6. We simulate 10^ paths for the estimation. 


Absolute error of the Monte Carlo MGF for model (V4) (log. scales) 



Figure 7.11. Absolute error (log-log plot) of the Monte Carlo estimate 
of the moment generating function E for model (Mdk defined in 


Section 7.5 The true values are estimated by the intercept of the linear 


extrapolation of the Monte Carlo estimates. The errors are calculated 
as the absolute difference between the intercept and the estimates. 10® 
paths were used in the simulation. 


the covariations observed in the market. For instance, the 33 x 33 covariation matrix of 
market yields with times to maturity ranging from 3 months to 30 years typically has 
rank between 5 and 9, as shown in|Figure 7.13 The ranks produced by the CRC models 
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o 


di 

o 


Absolute error of the Monte Carlo MGF for model (CIR4) (log. scales) 



Figure 7.12. Absolute error (log-log plot) of the Monte Carlo estimate 
of the moment generating function E for model (CIH4|. The 

true values are estimated by the intercept of the linear extrapolation of 
the Monte Carlo estimates. The errors are calculated as the absolute 
difference between the intercept and the estimates for different values 
of 6. 10^ paths were used in the simulation. 



(V1) 

(V2) 

(V3) 

(V4) 

(CIR1) 

(CIR2) 

(CIR3) 

{CIR4) 

Mean (%) 

0.25 

0.24 

0.25 

0.24 

0.25 

0.25 

0.25 

0.24 

Median (%) 

0.25 

0.24 

0.25 

0.24 

0.20 

0.12 

0.14 

0.19 

Volatility (%) 

0.09 

0.15 

0.13 

0.07 

0.19 

0.33 

0.30 

0.20 

Skewness 

-0.02 

0.00 

0.02 

0.00 

1.79 

2.75 

2.63 

1.90 

Kurtosis 

3.01 

2.99 

4.94 

3.21 

8.24 

14.47 

13.61 

9.06 

Minimum (%) 

-0.19 

-0.36 

-0.82 

-0.22 

0.00 

0.00 

0.00 

0.00 

Maximum (%) 

0.60 

0.82 

1.22 

0.61 

2.06 

4.61 

3.95 

2.21 

1st Ouartiie (%) 

0.18 

0.14 

0.17 

0.20 

0.11 

0.04 

0.05 

0.10 

3rd Quartile (%) 

0.31 

0.34 

0.32 

0.29 

0.33 

0.33 

0.33 

0.33 


Table 7. 1. Statistics of the short rate r(l) at time 1 in the models 
defined in Section 7.5 obtained by Monte-Carlo simulations. 10® paths 
and a step size of <5 = 0.02 were used in the simulation. 


(\Q and (CIlQ typically lie between 3 and 5. Thus, they are higher than those of the 
non-CRC models, but not as high as those of the market. 

Numerically, the covariation matrix is calculated as in (5.111 and the ranks are defined 
as the number of singular values differing significantly from zero. A comparison with 
Figures [7.15| and |7.14| shows that higher ranks are also related to higher volatility of the 
parameters. 
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Rank of covariation matrix 



Figure 7.13. Historical rank of the empirical covariation matrix (5.11) 
based on time windows oi M = 100 market yields with 33 different times 
to maturity n G {0.25, 0.5, 0.75,1,2,3,..., 30}. For comparison, the plot 
also features the average ranks obtained in simulations of the Hull-White 
extended affine models (\§ and (CIlQ as well as their CRC counterparts 
(Mi and (CIlQ. These models were calibrated using time windows of 
M = 100 observations. The missing values in (CIlQ and (CIlQ are due 
to non-admissible negative levels of mean reversion at these dates, see 
Section 7.4 and Figure 7.9 The averages are taken over 10^ simulated 


paths. In the numerical computation of the rank, eigenvalues which are 
10“® times smaller than the largest eigenvalue are rounded down to zero. 


Historical volatility of p 



time (t) 


Figure 7.14. Historical values of the parameter ai in the models (M^ 
and (CIlQ defined in [Section 7.5| estimated using time windows of 
M = 100 observations. 
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Historical volatiiity of a and a 



Figure 7.15. Historical values of the parameter <72 in the models (\Q 
and (CIlQ defined in Section 7.5 estimated using time windows of 
M = 100 observations. 


Appendix A. Consistent recalibration of Cox-Ingersoll-Ross models 

A.l. Overview. We describe CRC models based on CIR short rates, giving a detailed 
description of the simulation and calibration schemes. For comparison, we briefly digress 
to the CIR++ model and its CRC version. 


A.2. Setup and notation. We use the setup of Sections 3.2 and 4.2 setting X = M_|_, £ = 
0, A = 1. We do not specify the parameter space Y, yet, but we assume that for each {x, y) G 
X X Y, the volatility and drift coefficients are given by Ay{x) = ayX and By{x) = j5yX for 
some ay G (0, 00 ) and /3y G (— 00 , 0). For simplicity, we again choose equidistant grids of 
times tn = nS and times to maturity r„ = nS, for all n G No, where <5 is a positive constant. 


A.3. Hull-White extended Cox-Ingersoll-Ross models. For each fixed set of pa¬ 
rameters (y, 0) G Y X C(K+), the SDE for the short rate process is 


(A.l) 


dr{t) = [0{t) + Pyr{t))dt + ayr{t)dW{t), 


where W is one-dimensional (J^(t))t>o-Brownian motion. Thus, (y,0) satisfies Assump¬ 
tion 3.1 if and only if 0(t) > 0, for all t G K+. The functional characteristics (F, R) from 


Section 3.41 are 

Fy{u) = 0, Ry{u) = '^u^ + j3yU, for all u G M. 

Letting 7 j, = jdy + 2ay, the solutions of the corresponding Riccati equations are 

2(eT'A - 1) 


(A.2) 


>^y{t) = 0 , ^y{t) = 


ly Py _ 1 ) ’ 


for all t > 0. 
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Thus, by Theorem 3.3[ the forward rates in the Hull-White extended CIR model (A.ll 
with fixed parameters {y,0) are given by h{t) = 'Hy{S{t)9,r{t)), where 

x){t) = — f 0{s)"^'y{T — s)ds — (r)x, for all {x, r) G M x K_|_. 

Jo 

In contrast to the Vasicek model, the integral kernel is more complicated, 

2(e^«--l)7,e^«-(7,-/3,) 


Kir) = 


-2^ye'yy^ 


jy + 1 ) - /3y Ky'^ - 1 ) ( 7 ^ + 1 ) - /3y (e'^y'" - 1 )) 


2 > 


and there does not seem to be a closed-form expression for 6 = Cy{h,x). Instead, it must 
be calculated numerically as described in [Section 3.8[ 

The HJM drift and volatility are 


(A.3) 


HJM 


{x){t) = Kir)'^vir)ayX, 


<ry {x){T) = -JddyxK{r): 


and the HJM equation for forward rates reads as 

(A.4) dhit) = iAh{t) + fi^^^ih{t){0))) dt + a^^^{h{t){0))dW{t). 

A.4. Cox-Ingersoll-Ross CRC models. As the factor process is a function of the 
forward rate process (i.e., X{t) = r{t) = /i(t,0)), the corresponding CRC models can 
be characterised by the process {h,Y) instead of {h,X,Y). Thus, in accordance with 
Theorem 4.5 and Definition 4.6[ a process {h, Y) with values in H x Y may be called a 


CRC model if h satisfies the SPDE 
(A.5) dh{t) = ^Ah{t) + dt + aYf^{h{t){0))dW{t), 


with drift and volatility crp™ as defined in Equation (A.3) To ensure that the drift 
and volatility are well-defined, for all t G K+, it must be assumed that h{t){0) > 0 holds, 
for all t G IR+ . The maximally admissible set I in the CIR model is exactly characterised 
by this condition. 


A.5. Simulation of Cox-Ingersoll-Ross CRC models. The CRC model is simulated 
as described in [Algorithm 4.2 Discretisation in time and time to maturity is done as for 
the Vasicek model. However, in contrast to the Vasicek model, simulating the short rate 
process and calibrating Hull-White extensions is done by numerical approximations of 
order two. The resulting algorithm is presented below. 


Algorithm A.l (Simulation). Given an initial curve of forward rates h{0) and the 
parameter process V, execute iteratively the following steps, for each n G Nq: 


(i) The values of 0(t„) = Cy((^)(h(t„),r(t„)) at times to maturity 0 and J are calculated 

to g = -h{tn) - 


by applying 


Lemma 3.8 


6»(t„)(0) = Ah{tn)iO) - /3y(t„)/l(ti)(0), 

0(t„)(J) = ^(h(t„)(J) + VI/P(,„)(J) r(t„)) + e{tn)iO). 

(ii) A random draw r(t„_|_i) = ^(t„+i) of the CIR process with parameter 

Y{tn) and time-dependent drift 9{tn) is created using the second-order scheme of [^. 
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(iii) {h(tn+i), Ahltn+i)) is calculated from [h{tn), Ah{tn),r{tn+i)) using 
with integrals approximated by the trapezoid rule: 


Lemma 4.3 


/i(t„+i)(T) = h{tn)i6 + r) + 4'r(t„)(<5 + r) r(t„) - r(t„+i) 

+ 2 + r) + , 

Ah(t„+i){T) = Ah{tn){S + r) + + T)r{tn) - d'y((^)(T) r(t„+i) 

+ 2 + r) + 0(t„)((5)'I'y((_^)(r)^ . 


Here, h(tn+i) must be calculated at all times to maturity whereas Ah{tn+i) is 
needed only at Tq = 0. 


A. 6 . Calibration of Cox-Ingersoll-Ross CRC models. We proceed as in the Vasicek 
case described in Section SToj with [Equation (5.11) replaced by 


(A.6) 


[r{-,Ti),r{-,Tj)]{tn) - [r{-,Ti),r{-,Tj)]{tn-M) _ 4'(rj) 


The function ik depends on a, [3 as shown in Equation (A.2) 


A.7. CIRH—h models in the CRC framework. In the CIR++ model Section 3.9], 
also known as deterministic shift-extended CIR model, the short rate process is defined by 
r{t) = X{t)+0(t), where A is a CIR process and 0 is a deterministic function of time. Note 
that this is a different time-inhomogeneity than the one described in jSection A.3| In partic¬ 
ular, the factor process X is time-homogeneous and does not coincide with the short rate. 

Forward rate curves are given by 


h{t) = S(t)0 - by'l’y - %X{t), 


where '^y is the same as in the CIR case. 


see 


Equation (A.2) Given the parameter vector 


y and the factor A, this equation allows to calibrate 0 to a given yield curve without 
having to invert a Volterra integral operator. The HJM equation of the CIR-I--I- model is 

dh{t) = (A/i(t) + M™(A(t)))dt + a™(A(t))dW(t), 

(A.7) ! - 

dX{t) = {by -I- PyX{t)) dt ^ayX{t)dW{t), 


where and are the same as in the CIR case, see Equation (A.3) 

The CRC extension of the CIR-I-+ model is obtained by replacing the constant parame¬ 


ter vector y in (A.7) by a stochastic process {Y{t))t>o- The resulting equation is easier to 


handle than its CIR counterpart for two reasons. First, there are no boundary conditions 
on h. Indeed, 9 is allowed to assume negative values and can be calibrated to any forward 
rate curve. Thus, Equation (A.7) is defined on the entire space El x K+. Second, the SDE 
for A does not depend on h. Therefore, one can first solve for A, and then construct a 
mild solution h by stochastic convolution Section 6.1]: 


h{t)=S{t)h{0)+ / St.s9^\f){X{s))ds+ / 

Jo Jo 

The SDE for A is finite-dimensional. Therefore, existence and uniqueness of A can be 
shown by standard methods. For example, assuming that Y is independent of W, one can 


condition on Y and use results on time-inhomogeneous affine processes 18 to construct A 
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Simulation of the CRC model is analogue to [Algorithm 4.2[ The recalibration step is 
easier because no Volterra equation is involved. 

A disadvantage of the model is the presence of the hidden factor X. In contrast to the 
CIR version, X is not a function of the forward rate curve and cannot be directly observed. 
This is a challenge for calibration. We suggest an analogue approach to [Section 4.12 First, 
/3v(t)) ^Y{t}, and X(t) can be identified from the instantaneous covariation 


d[ri-,Ti),r{-,Tj)] (t) = aY{t) 




X(t)dt, 


of yields with times to maturity Ti,Tj. Subsequently, by^t) can be calibrated by least 
squares to the prevailing yield curve. Note that in this approach, X (t) is identified from the 
yield curve dynamics instead of extracted from the prevailing yield curve as in the Vasicek 
and CIR cases. For this reason the calibration is expected to be numerically more difficult. 
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